package ua.lviv.franko.solvers.stressanalysis.twoDimention.combined;

import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashSet;
import java.util.List;
import java.util.Set;

import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;
import org.jzy3d.maths.Range;
import org.jzy3d.plot3d.builder.Mapper;

import ua.lviv.franko.femShapeFunction2D.FEMShapeFunction1;
import ua.lviv.franko.integration.ConstantFunction2D;
import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.integration.VectorFunctionWithoutDerivatives2D;
import ua.lviv.franko.matrix.IRow;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.Solver;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.AbstractStressElement;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.LinearRectangularElement;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionLinear;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionStressInterface;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionType;
import ua.lviv.franko.triangulation.SaveStructures.Node;
import ua.lviv.franko.triangulation.SaveStructures.NodePair;
import ua.lviv.franko.triangulation.SaveStructures.Rectangle;
import ua.lviv.franko.visualisation.NewOkCancelDialog;
import ua.lviv.franko.visualisation.Visualise3D;

/**
 * @author marko
 * 
 *         this solver will use 2d models for both domains - main domain and slim inclusion
 * 
 */
public class CombinedSolverForMembrane2D extends Solver {

	protected static final double						EPS							= 0.000000000001;
	protected static final int							NUMBER_OF_UNKNOWNS			= 2;
	protected static final int							NUMBER_OF_QADRATURE_POINTS	= 9;

	protected List<AbstractStressElement>				elementsListSmall;
	protected List<BounderyConditionStressInterface>	bounderyConditionsSmall;
	protected HashMatrixLUStyle								bigMainMatrix;
	protected double[]									bigBVector;
	protected List<Node>								contactNodesList;
	protected Gauss2D									quadrature;
	protected List<AbstractStressElement>				nearBounderyElementsBig;
	protected List<AbstractStressElement>				nearBounderyElementsSmall;
	protected int										numBig;
	protected int										numSmall;

	public CombinedSolverForMembrane2D(List<Node> contactNodesList, int numBig, List<AbstractStressElement> elementsListBigProblem,
			List<BounderyConditionStressInterface> bounderyConditionsBigProblem, int numSmall, List<AbstractStressElement> elementsListSmallProblem,
			List<BounderyConditionStressInterface> bounderyConditionsSmallProblem) {
		super(numBig, elementsListBigProblem, bounderyConditionsBigProblem);
		this.numBig = numBig;
		this.numSmall = numSmall;
		this.elementsListSmall = elementsListSmallProblem;
		this.bounderyConditionsSmall = bounderyConditionsSmallProblem;
		this.contactNodesList = contactNodesList;
		quadrature = new Gauss2D(NUMBER_OF_QADRATURE_POINTS);
	}

	@Override
	public void generate() {
		super.generate(elementsList, bounderyConditions);
		super.applyBounderyConditions(bounderyConditions);
		// copy matrix for bigger problem;
		bigMainMatrix = mainMatrix;
		List<Double> bTmpVector = new ArrayList<Double>(b.length);
		for (int i = 0; i < b.length; i++) {
			bTmpVector.add(b[i]);
		}
		// copy for smaller problem
		final int shiftNumber = b.length;
		setN(numSmall);
		super.generate(elementsListSmall, bounderyConditionsSmall);
		super.applyBounderyConditions(bounderyConditionsSmall);
		List<IRow> rows = bigMainMatrix.getRows();
		List<IRow> tmpRows = mainMatrix.getRows();
		for (int i = 0; i < b.length; i++) {
			bTmpVector.add(b[i]);
			IRow row = tmpRows.get(i);
			List<Integer> pos = row.getIndexes();
			for (int j = 0; j < pos.size(); j++) {
				pos.set(j, pos.get(j) + shiftNumber);
			}
			rows.add(row);
		}

		// Lagrange constraints for displacements
		int rowsToAdd = 2;
		bigMainMatrix.addRows(rowsToAdd);
		for (int j = 0; j < rowsToAdd; j++) {
			bTmpVector.add(0.0);
		}

		for (int index = 0; index < contactNodesList.size() - 1; index++) {
			bigMainMatrix.addRows(rowsToAdd);
			for (int j = 0; j < rowsToAdd; j++) {
				bTmpVector.add(0.0);
			}
			Node A = contactNodesList.get(index);
			Node B = contactNodesList.get(index + 1);
			Node C = new Node(A.x + B.x, A.y + B.y, false, -1);
			final FEMShapeFunction1 fA = new FEMShapeFunction1(A, B, C);
			applyLagrange(0, 1.0, A, B, elementsList, fA);
			applyLagrange(shiftNumber, -1.0, A, B, elementsListSmall, fA);
		}

//		bigMainMatrix.addRows(6);
//		for (int i = 0; i < 6; i++) {
//			bTmpVector.add(0.0);
//		}
//
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 6, 13, 1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 5, 14, 1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 4, 15, 1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 3, 16, 1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 2, 17, 1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 1, 18, 1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 6, 19, -1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 5, 20, -1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 4, 21, -1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 3, 22, -1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 2, 23, -1.0);
//		bigMainMatrix.setValue(bigMainMatrix.getN() - 1, 24, -1.0);
//
//		bigMainMatrix.setValue(13, bigMainMatrix.getN() - 6, 1.0);
//		bigMainMatrix.setValue(14, bigMainMatrix.getN() - 5, 1.0);
//		bigMainMatrix.setValue(15, bigMainMatrix.getN() - 4, 1.0);
//		bigMainMatrix.setValue(16, bigMainMatrix.getN() - 3, 1.0);
//		bigMainMatrix.setValue(17, bigMainMatrix.getN() - 2, 1.0);
//		bigMainMatrix.setValue(18, bigMainMatrix.getN() - 1, 1.0);
//		bigMainMatrix.setValue(19, bigMainMatrix.getN() - 6, -1.0);
//		bigMainMatrix.setValue(20, bigMainMatrix.getN() - 5, -1.0);
//		bigMainMatrix.setValue(21, bigMainMatrix.getN() - 4, -1.0);
//		bigMainMatrix.setValue(22, bigMainMatrix.getN() - 3, -1.0);
//		bigMainMatrix.setValue(23, bigMainMatrix.getN() - 2, -1.0);
//		bigMainMatrix.setValue(24, bigMainMatrix.getN() - 1, -1.0);

		this.bigBVector = new double[bTmpVector.size()];
		for (int i = 0; i < bTmpVector.size(); i++) {
			this.bigBVector[i] = bTmpVector.get(i);
		}

//		for (int i = 0; i < bigMainMatrix.getN(); i++) {
//			for (int j = 0; j < bigMainMatrix.getN(); j++) {
//				System.out.print(bigMainMatrix.getValue(i, j) + "\t");
//			}
//			System.out.println();
//		}
//		System.out.println("======================");

	}

	private void applyLagrange(int shift, final double multiplier, Node A, Node B, List<AbstractStressElement> nearBounderyElements, final FEMShapeFunction1 fA) {
		for (final AbstractStressElement element : nearBounderyElements) {
			double[] results = quadrature.gaussOnGamma(A, B, new VectorFunctionWithoutDerivatives2D() {

				@Override
				public int vectorSize() {
					return element.getBaseFunctions().size() * 2;
				}

				@Override
				public double[] calculate(double x, double y) {
					double[] res = new double[vectorSize()];
					fA.setNodeNum(1);
					double faVal = fA.calculate(x, y) * multiplier;
					for (int i = 0; i < element.getBaseFunctions().size(); i++) {
						res[i] = faVal * element.getBaseFunctions().get(i).calculate(x, y);
					}
					fA.setNodeNum(2);
					faVal = fA.calculate(x, y) * multiplier;
					for (int i = 0; i < element.getBaseFunctions().size(); i++) {
						res[i] = faVal * element.getBaseFunctions().get(i).calculate(x, y);
					}
					return res;
				}
			});
			List<Integer> pos = element.getNodesNumbers();
			for (int i = 0; i < pos.size(); i++) {
				int row = bigMainMatrix.getN() - 4;
				int col = NUMBER_OF_UNKNOWNS * pos.get(i) + shift;
				bigMainMatrix.addValue(row, col, results[i]);
				bigMainMatrix.addValue(row + 1, col + 1, results[i]);

				bigMainMatrix.addValue(col, row, results[i]);
				bigMainMatrix.addValue(col + 1, row + 1, results[i]);

				row += 2;
				bigMainMatrix.addValue(row, col, results[i]);
				bigMainMatrix.addValue(row + 1, col + 1, results[i]);

				bigMainMatrix.addValue(col, row, results[i]);
				bigMainMatrix.addValue(col + 1, row + 1, results[i]);
			}
		}
	}

	public double[] solve() {
		LUFactorization lu = new LUFactorization(this.bigMainMatrix);
		return lu.calculate(bigBVector);
	}

	public static void main(String[] args) {

		double e = 0.6;
		double v = 0.6;
		IFunction2D f1 = new ConstantFunction2D(1.0);
		IFunction2D f3 = new ConstantFunction2D(1.0);

		final List<AbstractStressElement> elemsBig = new ArrayList<AbstractStressElement>();
		List<BounderyConditionStressInterface> bounderiesBig = new ArrayList<BounderyConditionStressInterface>();
		int n = 10;
		int m = 10;
		List<Node> nodes = new ArrayList<Node>();
		final double ax = 0.0, bx = 1.0, ay = 0.0, by = 1.0;
		prepare(ax, bx, ay, by, n, m, elemsBig, bounderiesBig, nodes, e, v, f1, f3);
		Set<Node> tmpNodes = new HashSet<Node>();
		for (AbstractStressElement element : elemsBig) {
			List<Integer> positions = element.getPos();
			Node A = nodes.get(positions.get(0));
			Node B = nodes.get(positions.get(1));
			Node C = nodes.get(positions.get(2));
			Node D = nodes.get(positions.get(3));
			// top
			if (Math.abs(A.y - by) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.y - by) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.y - by) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.y - by) < EPS) {
				tmpNodes.add(D);
			}
		}

		final List<AbstractStressElement> elemsSmall = new ArrayList<AbstractStressElement>();
		List<BounderyConditionStressInterface> bounderiesSmall = new ArrayList<BounderyConditionStressInterface>();
		int n2 = 10;
		int m2 = 10;
		final double ay2 = by;
		double by2 = by+0.05;
		double e2 = e;
		double v2 = v;
		prepare(ax, bx, ay2, by2, n2, m2, elemsSmall, bounderiesSmall, new ArrayList<Node>(), e2, v2, f1, f3);
		List<Node> contactNodesList = new ArrayList<Node>();
		contactNodesList.addAll(tmpNodes);
		Collections.sort(contactNodesList, new Comparator<Node>() {

			@Override
			public int compare(Node o1, Node o2) {
				return o1.x > o2.x ? 1 : -1;
			}
		});
		CombinedSolverForMembrane2D solver = new CombinedSolverForMembrane2D(contactNodesList, (n + 1) * (m + 1), elemsBig, bounderiesBig, (n2 + 1) * (m2 + 1), elemsSmall,
				bounderiesSmall);
		solver.generate();

		final double[] res = solver.solve();
		for (int j = 0; j < res.length; j++) {
			System.out.println(res[j]);
		}
		System.out.println();

//		Visualise3D.show(new Mapper() {
//			@Override
//			public double f(double x, double y) {
//				double val = 0.0;
//				int i = 0;
//				while (val == 0.0 && i < elemsBig.size()) {
//					val += elemsBig.get(i).calculateV(x, y, res);
//					i++;
//				}
//
//				// return ff.calculateDerivatyY(x, y);//(x, y);
//				return val;
//			}
//		}, new Range(ax, bx), new Range(ay, by));

		final double[] res2 = new double[(n2 + 1) * (m2 + 1) * NUMBER_OF_UNKNOWNS];
		for (int i = 0; i < res2.length; i++) {
			res2[i] = res[i + (n + 1) * (m + 1)*NUMBER_OF_UNKNOWNS];
		}
//		Visualise3D.show(new Mapper() {
//			@Override
//			public double f(double x, double y) {
//				double val = 0.0;
//				int i = 0;
//				while (val == 0.0 && i < elemsSmall.size()) {
//					val += elemsSmall.get(i).calculateV(x, y, res2);
//					i++;
//				}
//
//				// return ff.calculateDerivatyY(x, y);//(x, y);
//				return val;
//			}
//		}, new Range(ax, bx), new Range(ay2, by2));
//
		Visualise3D.show(new Mapper() {
			@Override
			public double f(double x, double y) {
				double val = 0.0;
				int i = 0;
				if (y < ay2) {
					while (val == 0.0 && i < elemsBig.size()) {
						val += elemsBig.get(i).calculateV(x, y, res);
						i++;
					}
				} else {
					while (val == 0.0 && i < elemsSmall.size()) {
						val += elemsSmall.get(i).calculateV(x, y, res2);
						i++;
					}
				}
				return val;
			}
		}, new Range(ax, bx), new Range(ay, by2));
		Visualise3D.show(new Mapper() {
			@Override
			public double f(double x, double y) {
				double val = 0.0;
				int i = 0;
				if (y < ay2) {
					while (val == 0.0 && i < elemsBig.size()) {
						val += elemsBig.get(i).calculateU(x, y, res);
						i++;
					}
				} else {
					while (val == 0.0 && i < elemsSmall.size()) {
						val += elemsSmall.get(i).calculateU(x, y, res2);
						i++;
					}
				}
				return val;
			}
		}, new Range(ax, bx), new Range(ay, by2));
		
		XYSeries s1 = new XYSeries("result1", true, true);
		XYSeries s2 = new XYSeries("result2", true, true);
		
		double myh = (bx-ax)/100.0;
		
		for (int i = 0; i < 100; i++) {
			double x = ax + i*myh;
			double y = 1.0;
			double val1 = 0.0;
			double val2 = 0.0;
			int j = 0;
			while (val1 == 0.0 && j < elemsSmall.size()) {
				val1 += elemsSmall.get(j).calculateU(x, y, res2);
				j++;
			}
			j = 0;
			while (val2 == 0.0 && j < elemsSmall.size()) {
				val1 += elemsSmall.get(j).calculateV(x, y, res2);
				j++;
			}
			s1.add(x, val1);
			s2.add(x, val2);
		}


		XYSeriesCollection col = new XYSeriesCollection();
		col.addSeries(s1);
		col.addSeries(s2);
		NewOkCancelDialog d = new NewOkCancelDialog(null, false, col);
		d.setLocationRelativeTo(null);
		d.setVisible(true);
	}

	public static void prepare(double ax, double bx, double ay, double by, int n, int m, List<AbstractStressElement> elems, List<BounderyConditionStressInterface> bounderies,
			List<Node> nodes, double e, double v,IFunction2D f1, IFunction2D f3) {

		double hx = (bx - ax) / (double) n, hy = (by - ay) / (double) m;

		Set<Node> tmpNodesSet = new HashSet<Node>();
		double startingAy = ay;
		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n; j++) {
				Rectangle r = new Rectangle();
				r.A = new Node(ax + j * hx, ay, false, i * (n + 1) + j);
				r.B = new Node(ax + (j + 1) * hx, ay, false, i * (n + 1) + j + 1);
				r.C = new Node(ax + (j + 1) * hx, ay + hy, false, (i + 1) * (n + 1) + j + 1);
				r.D = new Node(ax + j * hx, ay + hy, false, (i + 1) * (n + 1) + j);
				tmpNodesSet.add(r.A);
				tmpNodesSet.add(r.B);
				tmpNodesSet.add(r.C);
				tmpNodesSet.add(r.D);

				LinearRectangularElement elem = new LinearRectangularElement(r.A, r.B, r.C, r.D, e, v, f1, f3);
				elems.add(elem);
			}
			ay += hy;
		}
		ay = startingAy;
		nodes.addAll(tmpNodesSet);
		Collections.sort(nodes);
		List<NodePair> leftNodePairs = new ArrayList<NodePair>();
		List<NodePair> rightNodePairs = new ArrayList<NodePair>();
		List<NodePair> topNodePairs = new ArrayList<NodePair>();
		List<NodePair> bottomNodePairs = new ArrayList<NodePair>();
		for (AbstractStressElement element : elems) {
			List<Integer> positions = element.getPos();
			Node A = nodes.get(positions.get(0));
			Node B = nodes.get(positions.get(1));
			Node C = nodes.get(positions.get(2));
			Node D = nodes.get(positions.get(3));
			// left
			List<Node> tmpNodes = new ArrayList<Node>();
			if (Math.abs(A.x - ax) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.x - ax) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.x - ax) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.x - ax) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				leftNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// right
			tmpNodes.clear();
			if (Math.abs(A.x - bx) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.x - bx) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.x - bx) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.x - bx) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				rightNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// top
			tmpNodes.clear();
			if (Math.abs(A.y - by) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.y - by) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.y - by) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.y - by) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				topNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// bottom
			tmpNodes.clear();
			if (Math.abs(A.y - ay) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.y - ay) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.y - ay) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.y - ay) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				bottomNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}
		}
		Gauss2D gauss2d = new Gauss2D(7);

		bounderies.add(new BounderyConditionLinear(BounderyConditionType.DIRIHLET_TYPE, new ConstantFunction2D(0.0), new ConstantFunction2D(0.0), leftNodePairs, gauss2d));

		bounderies.add(new BounderyConditionLinear(BounderyConditionType.DIRIHLET_TYPE, new ConstantFunction2D(0.0), new ConstantFunction2D(0.0), rightNodePairs, gauss2d));

	}
}